library(xtable)
library("arm")
library("foreign")
library("blme")
library("texreg")
library("survey")

setwd(".../Final replication file")


# replicating the original artcile's results based on code from Warshaw and Rodden (2012)
source("source_MRP_US.R")

# District estimates and multilevel model:
replication.model <- marriage.model



#########################################################################
#### Modified Replication Code from Warshaw and Rodden (2012) ###########
#########################################################################
#Analysis for State Senate Districts

# District id variable
districts<-read.csv(paste("districts_ssd.csv",sep=""))
district.uid<-districts$ssd_district_uid


### Prepare survey data for response model with state senate variables 
data1<-read.dta("gay_marriage_senate.dta")
state.gaymarriage<-data1$state_gaymarriage
race2<-data1$race
cst2<-factor(data1$cst)
education5<-factor(data1$education)
factor.education<-factor(data1$education5)
region<-data1$region_num
districturban.full<-data1$ssd_district_urban
districtincome.full<-data1$ssd_district_income/100000
districtveteran.full<-data1$ssd_district_veteran
statereligion.full<-data1$state_religion
stateunion.full<-data1$state_union
districtgay.full<-data1$ssd_district_gay
district_num<-data1$ssd_district_uid
gender<-data1$gender

# adding age as predictor
age <- data1$age_exact
age.group <- rep(0,length(age))
age.group[age<20] <- 0
age.group[age>19&age<25] <- 1
age.group[age>24&age<30] <- 2
age.group[age>29&age<35] <- 3
age.group[age>34&age<40] <- 4
age.group[age>39&age<45] <- 5
age.group[age>44&age<50] <- 6
age.group[age>49&age<55] <- 7
age.group[age>54&age<60] <- 8
age.group[age>59&age<65] <- 9
age.group[age>64&age<70] <- 10
age.group[age>69&age<75] <- 11
age.group[age>74&age<80] <- 12
age.group[age>79&age<85] <- 13
age.group[age>84&age<90] <- 14
age.group[age>89] <- 15



##############################################################################################################
##### Post-Stratification Step
new.data <- read.dta("DATA_v4_ageshares.dta")
new.data1 <- new.data[,c(47,53:67)] 
colnames(new.data1)

## Survey data

DATAmodel <- data.frame(cbind(age.group, state.gaymarriage, race2, gender, cst2, district_num, education5, districtincome.full, districturban.full, districtveteran.full, statereligion.full, stateunion.full, region, districtgay.full))


# generating sample information:
# race (2-3-1-4)- gender (0-1) - education (1:5) - AGE.GROUP (1:15)
table(DATAmodel$race2)
table(DATAmodel$gender)
table(DATAmodel$education5)
DATAmodel$age.group
AGE <- factor(DATAmodel$age.group, levels = c(0:15))
Survey.Age.Info <- matrix(NA,40,16)
Survey.Age.Info[1,] <- table(AGE[DATAmodel$race2==2 & DATAmodel$gender==0 & DATAmodel$education5==1])
Survey.Age.Info[2,] <- table(AGE[DATAmodel$race2==2 & DATAmodel$gender==0 & DATAmodel$education5==2])
Survey.Age.Info[3,] <- table(AGE[DATAmodel$race2==2 & DATAmodel$gender==0 & DATAmodel$education5==3])
Survey.Age.Info[4,] <- table(AGE[DATAmodel$race2==2 & DATAmodel$gender==0 & DATAmodel$education5==4])
Survey.Age.Info[5,] <- table(AGE[DATAmodel$race2==2 & DATAmodel$gender==0 & DATAmodel$education5==5])
Survey.Age.Info[6,] <- table(AGE[DATAmodel$race2==2 & DATAmodel$gender==1 & DATAmodel$education5==1])
Survey.Age.Info[7,] <- table(AGE[DATAmodel$race2==2 & DATAmodel$gender==1 & DATAmodel$education5==2])
Survey.Age.Info[8,] <- table(AGE[DATAmodel$race2==2 & DATAmodel$gender==1 & DATAmodel$education5==3])
Survey.Age.Info[9,] <- table(AGE[DATAmodel$race2==2 & DATAmodel$gender==1 & DATAmodel$education5==4])
Survey.Age.Info[10,] <- table(AGE[DATAmodel$race2==2 & DATAmodel$gender==1 & DATAmodel$education5==5])
Survey.Age.Info[11,] <- table(AGE[DATAmodel$race2==3 & DATAmodel$gender==0 & DATAmodel$education5==1])
Survey.Age.Info[12,] <- table(AGE[DATAmodel$race2==3 & DATAmodel$gender==0 & DATAmodel$education5==2])
Survey.Age.Info[13,] <- table(AGE[DATAmodel$race2==3 & DATAmodel$gender==0 & DATAmodel$education5==3])
Survey.Age.Info[14,] <- table(AGE[DATAmodel$race2==3 & DATAmodel$gender==0 & DATAmodel$education5==4])
Survey.Age.Info[15,] <- table(AGE[DATAmodel$race2==3 & DATAmodel$gender==0 & DATAmodel$education5==5])
Survey.Age.Info[16,] <- table(AGE[DATAmodel$race2==3 & DATAmodel$gender==1 & DATAmodel$education5==1])
Survey.Age.Info[17,] <- table(AGE[DATAmodel$race2==3 & DATAmodel$gender==1 & DATAmodel$education5==2])
Survey.Age.Info[18,] <- table(AGE[DATAmodel$race2==3 & DATAmodel$gender==1 & DATAmodel$education5==3])
Survey.Age.Info[19,] <- table(AGE[DATAmodel$race2==3 & DATAmodel$gender==1 & DATAmodel$education5==4])
Survey.Age.Info[20,] <- table(AGE[DATAmodel$race2==3 & DATAmodel$gender==1 & DATAmodel$education5==5])
Survey.Age.Info[21,] <- table(AGE[DATAmodel$race2==1 & DATAmodel$gender==0 & DATAmodel$education5==1])
Survey.Age.Info[22,] <- table(AGE[DATAmodel$race2==1 & DATAmodel$gender==0 & DATAmodel$education5==2])
Survey.Age.Info[23,] <- table(AGE[DATAmodel$race2==1 & DATAmodel$gender==0 & DATAmodel$education5==3])
Survey.Age.Info[24,] <- table(AGE[DATAmodel$race2==1 & DATAmodel$gender==0 & DATAmodel$education5==4])
Survey.Age.Info[25,] <- table(AGE[DATAmodel$race2==1 & DATAmodel$gender==0 & DATAmodel$education5==5])
Survey.Age.Info[26,] <- table(AGE[DATAmodel$race2==1 & DATAmodel$gender==1 & DATAmodel$education5==1])
Survey.Age.Info[27,] <- table(AGE[DATAmodel$race2==1 & DATAmodel$gender==1 & DATAmodel$education5==2])
Survey.Age.Info[28,] <- table(AGE[DATAmodel$race2==1 & DATAmodel$gender==1 & DATAmodel$education5==3])
Survey.Age.Info[29,] <- table(AGE[DATAmodel$race2==1 & DATAmodel$gender==1 & DATAmodel$education5==4])
Survey.Age.Info[30,] <- table(AGE[DATAmodel$race2==1 & DATAmodel$gender==1 & DATAmodel$education5==5])
Survey.Age.Info[31,] <- table(AGE[DATAmodel$race2==4 & DATAmodel$gender==0 & DATAmodel$education5==1])
Survey.Age.Info[32,] <- table(AGE[DATAmodel$race2==4 & DATAmodel$gender==0 & DATAmodel$education5==2])
Survey.Age.Info[33,] <- table(AGE[DATAmodel$race2==4 & DATAmodel$gender==0 & DATAmodel$education5==3])
Survey.Age.Info[34,] <- table(AGE[DATAmodel$race2==4 & DATAmodel$gender==0 & DATAmodel$education5==4])
Survey.Age.Info[35,] <- table(AGE[DATAmodel$race2==4 & DATAmodel$gender==0 & DATAmodel$education5==5])
Survey.Age.Info[36,] <- table(AGE[DATAmodel$race2==4 & DATAmodel$gender==1 & DATAmodel$education5==1])
Survey.Age.Info[37,] <- table(AGE[DATAmodel$race2==4 & DATAmodel$gender==1 & DATAmodel$education5==2])
Survey.Age.Info[38,] <- table(AGE[DATAmodel$race2==4 & DATAmodel$gender==1 & DATAmodel$education5==3])
Survey.Age.Info[39,] <- table(AGE[DATAmodel$race2==4 & DATAmodel$gender==1 & DATAmodel$education5==4])
Survey.Age.Info[40,] <- table(AGE[DATAmodel$race2==4 & DATAmodel$gender==1 & DATAmodel$education5==5])


colnames(Survey.Age.Info) <- c("<20", "20--24","25--29", "30--34", "35--39", "40--44", "45--49", "50--54", "55--59", "60--64","65--69","70--74","75--79","80--84","85--90",">90")
print(xtable(Survey.Age.Info),include.rownames=FALSE)
# We will later need to know how AGE was distributed in the survey data
Survey.Age.frequency <- colSums(Survey.Age.Info[,-1])/sum(Survey.Age.Info[,-1])

tabL <- matrix(NA,16,3)
tabL[2:16,1] <- round(Survey.Age.frequency,3)
tabL[2:16,2] <- unlist(round(new.data1[1,-1],3))
tabL[2:16,3] <- round(tabL[2:16,2]/tabL[2:16,1],3)

xtable(tabL)


# Read in census data without age
### Prepare data for predictions with state senate data (n=77680, which is 2*4*5 * 1942) 
PS.Data<-read.csv(paste("Post-Stratification Data - RImport (ssd).csv",sep=""))



# CHANGE TO NEW SET-UP; 2 * 4 * 5 * 15 * 1942
# race - gender - education - age
b <- rep(15,dim(PS.Data)[1])
PS.Data1 <- PS.Data[rep(1:nrow(PS.Data), b), ]



for (a in 1:1942){
		a1 <- (a-1)*600 + 1	
		a2 <- a*600
		district.census.Age.frequency <- new.data1[a,-1]
		# create ratio of actual age group distribution and survey distribution
		district.corrector <- district.census.Age.frequency/Survey.Age.frequency
		Age.district <- Survey.Age.Info[,-1] %*% diag(district.corrector)
		Age.district.percent <- Age.district/rowSums(Age.district)
		Age.district.percent.vec <- c(t(Age.district.percent))
		# this creates the relative proportions per cell via "adjusted synthetic" distribution
		PS.Data1[c(a1:a2),6] <- PS.Data1[c(a1:a2),6] * Age.district.percent.vec
		print(a)
}


# now, that I adapted the proportions, I need to add the age category
age.category <- rep(c(1:15),77680)
PS.Data1 <- cbind(PS.Data1, age.category)

proportion<-PS.Data1$proportion
ps_cregion<-PS.Data1$ps_cregion
cst_ps<-PS.Data1$cst_ps
cat_educ<-PS.Data1$cat_educ
cat_race<-PS.Data1$cat_race
cat_female<-PS.Data1$cat_female
ps_district_income <-PS.Data1$ps_district_income/100
ps_district_urban<-PS.Data1$ps_district_urban
ps_district_veteran<-PS.Data1$ps_district_veteran
ps_districtgay<-PS.Data1$ps_district_gay
ps_state_religion<-PS.Data1$ps_state_religion
ps_state_union<-PS.Data1$ps_state_union
ps_district_uid<-PS.Data1$ssd_district_uid
cat_age<-PS.Data1$age.category


###########	
# Prepare data for disaggregation (calculate mean for each senate district)
sample.mean.gaymarriage<-tapply(state.gaymarriage, data1$ssd_district_uid, mean)
sample.count.gaymarriage<-tapply(state.gaymarriage, data1$ssd_district_uid, length)
sample.mean.gaymarriage.names<-as.numeric(rownames(sample.mean.gaymarriage))
sample.mean.gaymarriage.adj1<-cbind(sample.mean.gaymarriage, sample.mean.gaymarriage.names)
sample.mean.gaymarriage.adj2<-data.frame(matrix(NA, nrow=1942, ncol=4))

sample.count.gaymarriage.names<-as.numeric(rownames(sample.count.gaymarriage))
sample.count.gaymarriage.adj1<-cbind(sample.count.gaymarriage, sample.count.gaymarriage.names)
# 1 senate district is missing


## Create a vector of district effects without any missing cells (zeros where no data for a district)
for (i in 1:length(unique(district.uid)))
{
sample.mean.gaymarriage.adj2[i,1]<-sum(sample.mean.gaymarriage.adj1[sample.mean.gaymarriage.names==i,1], na.rm=FALSE)
sample.mean.gaymarriage.adj2[i,2]<-sum(sample.count.gaymarriage.adj1[sample.count.gaymarriage.names==i,1], na.rm=FALSE)
}

sample.mean.gaymarriage.adj2[,3]<-districts$cst

sample.mean.gaymarriage.adj2[,4]<-sample.mean.gaymarriage.adj2[,1]

for (i in 1:length(unique(district.uid))) {
				if (sample.mean.gaymarriage.adj2[i,2]==0){
									sample.mean.gaymarriage.adj2[i,4]<-NA
									}
				}
sample.mean.gaymarriage<-sample.mean.gaymarriage.adj2[,4]



#########################################################################################################
## Marriage response model with AGE as explanatory variable (w/ AGE)
#########################################################################################################

marriage.model.wAGE <-  glmer(formula = state.gaymarriage ~ (1 | age.group) + (1 | race2) + gender + (1 | cst2) + (1 | district_num) +  (1 | education5)+ districtincome.full + districturban.full + districtveteran.full + statereligion.full + stateunion.full + (1 | region) + districtgay.full, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa"))


# Get RE for each district (1 is missing)
district.effects<-ranef(marriage.model.wAGE)$district_num
district.effects.names<-as.numeric(rownames(district.effects))
district.effects2<-cbind(district.effects, district.effects.names)

district.effects3<-data.frame(matrix(NA, nrow=1942, ncol=1))

# Create a vectors of district effects without any missing cells (zeros where no data for a district)
for (i in 1:length(unique(district.uid))){
	j<-"i"
	district.effects3[i,1]<-sum(district.effects2[district.effects2$district.effects.names==i,1], na.rm=TRUE)
}


# Get RE for each state 
state.effects<-ranef(marriage.model.wAGE)$cst2
state.effects.names<-rownames(state.effects)
state.effects2<-cbind(state.effects, state.effects.names)
state.effects3<-data.frame(matrix(NA, nrow=51, ncol=2))
states<-unique(cst_ps)
state.effects3[,2]<-states

state.effects3[1:48,1]<-state.effects2[,1]
states2<-factor(states)

## Create a vector of district effects without any missing cells (zeros where no data for a district)
for (i in 1:length(states)){
	state.effects3[i,1]<-sum(state.effects2[state.effects.names==states[i],1], na.rm=TRUE)
}
state.effects3[50,2]<-"AK"
state.effects3[51,2]<-"HI"

# Make predictions
cell.pred<-invlogit(fixef(marriage.model.wAGE)["(Intercept)"]+ ranef(marriage.model.wAGE)$age.group[cat_age,1] +ranef(marriage.model.wAGE)$education5[cat_educ,1] + ranef(marriage.model.wAGE)$race2[cat_race,1] +   fixef(marriage.model.wAGE)["gender"] *cat_female +    district.effects3[ps_district_uid,1] + state.effects3[cst_ps,1]  +  ranef(marriage.model.wAGE)$region[ps_cregion,1]  + fixef(marriage.model.wAGE)["districtincome.full"]*ps_district_income + fixef(marriage.model.wAGE)["districturban.full"]*ps_district_urban + fixef(marriage.model.wAGE)["districtveteran.full"]*ps_district_veteran+fixef(marriage.model.wAGE)["statereligion.full"]*ps_state_religion + fixef(marriage.model.wAGE)["stateunion.full"]*ps_state_union+ fixef(marriage.model.wAGE)["districtgay.full"]*ps_districtgay)




# Post-stratification
cell.pred.weighted <- cell.pred * proportion

mrp.mean.gaymarriage<-tapply(cell.pred.weighted, PS.Data1$ssd_district_uid, sum)
write.csv(mrp.mean.gaymarriage,file="stategaymarriage.csv")

mrp.mean.gaymarriage2<-data.frame(matrix(NA, nrow=1942, ncol=3))
mrp.mean.gaymarriage2<-cbind(mrp.mean.gaymarriage, districts$ssd_fips_num)



# Get disaggregation and MrP data for the congressional districts for each of the 5 states as well as the referendum results from the datasets

ohio <-read.csv(paste("OH_GM_SSD.csv",sep=""))
mrp.ohio<-mrp.mean.gaymarriage2[1389:1421,]
raw.ohio<-data.frame(matrix(NA, nrow=33, ncol=1))
raw.ohio<-sample.mean.gaymarriage[1389:1421]

ca <-read.csv(paste("CA_GM_SSD.csv",sep=""))
mrp.ca<-mrp.mean.gaymarriage2[121:160,]
raw.ca<-data.frame(matrix(NA, nrow=33, ncol=1))
raw.ca<-sample.mean.gaymarriage[121:160]

wi<-read.csv(paste("WI_GM_SSD.csv",sep=""))
mrp.wi<-mrp.mean.gaymarriage2[1880:1912,]
raw.wi<-data.frame(matrix(NA, nrow=33, ncol=1))
raw.wi<-sample.mean.gaymarriage[1880:1912]

mi <-read.csv(paste("MI_GM_SSD.csv",sep=""))
mrp.mi<-mrp.mean.gaymarriage2[815:852,]
raw.mi<-data.frame(matrix(NA, nrow=33, ncol=1))
raw.mi<-sample.mean.gaymarriage[815:852]

az <-read.csv(paste("AZ_GM_SSD.csv",sep=""))
mrp.az<-mrp.mean.gaymarriage2[56:85,]
raw.az<-data.frame(matrix(NA, nrow=33, ncol=1))
raw.az<-sample.mean.gaymarriage[56:85]
raw.az[15]<-NA


ssd.disa<-append(raw.ohio, raw.ca)
ssd.disa<-append(ssd.disa, raw.wi)
ssd.disa<-append(ssd.disa, raw.mi)
ssd.disa<-append(ssd.disa, raw.az)
ssd.mrsp<-append(mrp.ohio[,1], mrp.ca[,1])
ssd.mrsp<-append(ssd.mrsp, mrp.wi[,1])
ssd.mrsp<-append(ssd.mrsp, mrp.mi[,1])
ssd.mrsp<-append(ssd.mrsp, mrp.az[,1])
ssd.referendum<-append(ohio[,2], ca[,2])
ssd.referendum<-append(ssd.referendum, wi[,2])
ssd.referendum<-append(ssd.referendum, mi[,2])
ssd.referendum<-append(ssd.referendum, az[,2])



# Squared error
se.dis <- (ssd.disa - ssd.referendum)^2						# disaggregation
se.mrsp <- (ssd.mrsp - ssd.referendum)^2					# MrsP approach
se.mrp.rep <- (ssd.mrp - ssd.referendum)^2		# MRP replication (Warshaw and Rodden, 2012)

# Mean squarred error
mse.dis <- mean(se.dis, na.rm=T)		# disaggregation
mse.mrsp <- mean(se.mrsp)				# MrsP approach
mse.mrp.rep <- mean(se.mrp.rep)			# MRP replication (Warshaw and Rodden, 2012)





# Plot

png("Figure04.png", width = 1800, height = 1000, pointsize=24)
#par(mfrow=c(1,1), family="CMU Serif") 
par(mfrow=c(1,1))
plot(ssd.mrsp, ssd.referendum, main="", xlab="", ylab="", axes=F, ylim=c(0,1),  xlim=c(0,1), pch=20, col="blue")
segments(0, 1,1,1,col="black", lty=1 )
segments(1, 0,1,1,col="black", lty=1 )
axis(1, at = c(0, .2, .4,.6,.8, 1),cex.axis=1.4, label=c("0%","20%", "40%",  "60%",  "80%",  "100%"), pos=0)
axis(2, at = c(0, .2, .4,.6,.8, 1),cex.axis=1.4, label=c("0%","20%", "40%",  "60%",  "80%",  "100%"), las = 1, tick = T, pos=0)
title(main="",xlab = "Predictions", ylab = "Vote Outcome", cex.lab=1.5)
clip(0, 1, 0, 1)
segments(0, 0,1,1,col="grey", lty=2 )
points(ssd.mrp, ssd.referendum, pch=2, col="red", cex=0.8, lwd=1.4)
points(ssd.disa, ssd.referendum, pch=3,col="grey40", cex=0.5)
legend(0.01,0.99,legend=c("MrsP", "MrP","Disaggregation"), col=c("blue","red","grey40"), pch=c(20,2,3), bty="n", pt.cex=c(1.5,.8,0.5),pt.lwd=c(NA,1.5,NA))

legend(0.63,0.18,legend=c(paste("MrsP MSE:", round(mse.mrsp,3)),paste("MrP MSE:", round(mse.mrp.rep,3)), paste("Disaggregation MSE", round(mse.dis,3))), col=c("blue","red","grey40"), pch=c(20,2,3), bty="n", pt.cex=c(1.5,.8,0.5),pt.lwd=c(NA,1.5,NA),cex=1.2, text.col=c("blue","red","grey40"))
dev.off()

#########################################################################################################
### For Appendix


# Plot of RE

test.sim <- sim(marriage.model.wAGE, n.sims=1000)


yy <- c(attributes(test.sim)$ranef$age.group)
xx <- rep(NA,length(yy))
for (i in 0:15){
	a <- 1000*i + 1
	b <- a +999
	xx[a:b] <- rep(i+1,1000)
}

png("Figure06.png", width = 1800, height = 1000, pointsize=24)
#par(mfrow=c(1,1), family="CMU Serif") 
par(mfrow=c(1,1))
plot(xx,yy, pch=19, col=rgb(204,204,204,04,maxColorValue=255), ylab="Random Effects for Age Groups",xaxt="n", xlab="Age Groups")
points(c(1:16), colMeans(attributes(test.sim)$ranef$age.group), pch=20)
axis(side=1,at=c(1:16), labels=c("<20", "20-24", "25-29", "30-34", "35-39", "40-44", "45-49", "50-54", "55-59", "60-64", "65-69", "70-74", "75-79", "80-84", "85-89", ">89"))
dev.off()

# Correlation among Level 1 Precitors

data.c <- data.frame(cbind(age, c(education5), race2, gender))
colnames(data.c) <- c("age", "education", "race",  "gender")

test.mat <- attributes(marriage.model.wAGE)$"frame"
test.mat <- cbind(test.mat$age,test.mat$education,test.mat$race,test.mat$gender)

cor(test.mat,use="pairwise.complete.obs")
xtable(cor(test.mat,use="pairwise.complete.obs"),caption = "Table 6 in Appendix")


# Output of Response Model (Warshaw and Rodden)


texreg(list(marriage.model.wAGE, replication.model), digits=2,  naive=TRUE, caption = "Table 7 in Appendix")






######################################################################################################
### Raking Rather Than Synthetic Post-Stratification
library("survey")


# age.shares (census info)
dim(new.data1)

### Prepare data for predictions with state senate data (n=77680, which is 2*4*5 * 1942) 
PS.Data<-read.csv(paste("Post-Stratification Data - RImport (ssd).csv",sep=""))
# data has proportions on 40 types per district and we will extend this to 600 types (to add 15 age categories)



# create fake indicator for SexGenderEduc-Variable (ranges from 1 to 40)
cat_RGE <- rep(c(1:40),15)

# add the indicator
SuperData <- cbind(PS.Data1,cat_RGE,cell.pred)


## We will now loop through each district and create a fake dataset with 600 voters. Then we use raking to 
# generate weights, then we look at the weighted average support.

mrr.estimates <- rep(NA, 1942)
for (i in 1:1942){
	d <- i
	d1 <- (d-1)*40 + 1	
	d2 <- d*40
	e1 <- (d-1)*600 + 1
	e2 <- d*600
	
	Super.Data.district.d <- SuperData[e1:e2,30:32]

	#This is district data from census
	marginal.RaceGenderEduc <- PS.Data[1:40,6]
	RGE.dist <- data.frame(cat_RGE = c(1:40), Freq = 600 * marginal.RaceGenderEduc)
	marginal.age <- unlist(new.data1[d,-1])
	age.dist <- data.frame(age.category = c(1:15), Freq = 600 * marginal.age)
	
	DATA.unweighted <- svydesign(ids=~1, data=Super.Data.district.d)
	DATA.svy.rake <- rake(design = DATA.unweighted, sample.margins = list(~age.category,~cat_RGE), population.margins = list(age.dist,RGE.dist))
	# District MR-Raking estimate
	mrr.estimates[i] <- svymean(~cell.pred,design=DATA.svy.rake)
}

summary(weights(DATA.svy.rake))

se.mrsp <- (ssd.mrsp - ssd.referendum)^2
se.mrr <- (mrr.estimates[c(1389:1421,121:160,1880:1912,815:852,56:85)] - ssd.referendum)^2

mse.mrsp <- mean(se.mrsp, na.rm=T)
mse.mrr <- mean(se.mrr, na.rm=T)

# Plot mit raking

png("Figure08.png", width = 1800, height = 1000, pointsize=24)
#par(mfrow=c(1,1), family="CMU Serif")
par(mfrow=c(1,1))
plot(ssd.mrsp, ssd.referendum, main="", xlab="", ylab="", axes=F, ylim=c(0,1),  xlim=c(0,1), pch=20, col=rgb(0,0,255,150,maxColorValue=255))
segments(0, 1,1,1,col="black", lty=1 )
segments(1, 0,1,1,col="black", lty=1 )
axis(1, at = c(0, .2, .4,.6,.8, 1),cex.axis=1.4, label=c("0%","20%", "40%",  "60%",  "80%",  "100%"), pos=0)
axis(2, at = c(0, .2, .4,.6,.8, 1),cex.axis=1.4, label=c("0%","20%", "40%",  "60%",  "80%",  "100%"), las = 1, tick = T, pos=0)
title(main="",xlab = "Predictions", ylab = "Vote Outcome", cex.lab=1.5)
clip(0, 1, 0, 1)
segments(0, 0,1,1,col="grey", lty=2 )
points(ssd.mrp, ssd.referendum, pch=20, col=rgb(255,0,0,150,maxColorValue=255))
points(mrr.estimates[c(1389:1421,121:160,1880:1912,815:852,56:85)], ssd.referendum, pch=20,col=rgb(160,32,240,150,maxColorValue=255), cex=1.0)
legend(0.01,0.99,legend=c("MrsP","MrR", "Classic MrP"), col=c("blue","purple","red"), pch=c(20,20,20), bty="n", pt.cex=c(1.5,1.5,1.5),cex=1.2)

	text(0.63,0.12,paste("MrsP MSE:", round(mse.mrsp,3)), pos=4, col="blue", cex=1.3)
	text(0.63,0.08,paste("MrR MSE: 0.010"), pos=4, col="purple", cex=1.3)
	text(0.63,0.04,paste("MrP MSE:", round(mse.mrp.rep,3)), pos=4, col="red", cex=1.3)	
dev.off()



 


